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Abstract. Recently, we have introduced a novel approach to deal with the sign problem that 
prevents the Monte Carlo simulations of a class of quantum field theories (QFTs). The idea is 
to formulate the QFT on a Lefschetz thimble. Here we review the formulation of our approach 
and describe the Aurora Monte Carlo algorithm that we are currently testing on a scalar field 
theory with a sign problem. 



1. Introduction 

Recently, we have proposed a novel approach [Tj to deal with the sign problem that hinders 
Monte Carlo simulations of many quantum field theories (QFTs). The approach consists in 
formulating the QFT on a Lefschetz thimble. In this paper we concentrate on the application 
to a scalar field theory with a sign problem. In particular, we review the formulation and the 
justification of the approach, and we also describe the Aurora Monte Carlo algorithm that we 
are currently testing. 



2. Formulation 

In this paper we consider a scalar field theory with U(l) global symmetry: 
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where 4>(x) is a complex scalar field. This model can be regularized on a lattice as: 
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and has a sign problem when [i ^ 0, which has been successfully treated in [2] and also in [3]. 
Our approach consists in defining the observables as: 
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where the integration domain J§ is the Lefschetz thimble [H [3] attached to a global minimum 
^giob °f the real part of the action Sr = KS, when restricted to the original domain C = M. 2V . 
More precisely, Jq is the manifold of real dimension n = 2V, defined as union of all the curves 
of steepest descent (SD) for Sr, i.e., solutions of: 
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that end in </> g i b for r — >• oo. Here we have assumed the usual complexification of the complex 
scalar field (see, e.g., [3]), where both the real (</>i) and imaginary (^2) part of the original fields 
become complex 4> a = <p£® + i<j)a\ a = 1,2. 



2.1. Formulation in presence of SSB 

Note that the formulation described above is possible only when c/> g i b is a minimum of the action 
with positive definite Hessian. In presence of symmetries that act non-trivially in <^>giob this is 
not the case. 

Gauge symmetries can be treated as described in [HJ EJ [T], where, essentially, the thimble 
is defined modulo gauge transformations. But this is not suitable to study the possibility of 
spontaneous symmetry breaking (SSB). However, the proper way to study SSB (even in the 
ordinary formulation of the functional integral) is to introduce a small term of explicit symmetry 
breaking and study the limit when such term goes to zero. In this way, 4> g i b becomes a true 
minimum and the thimble can be defined. 



2.2. Justification of the approach 

The functional integral in Eq. ^ does not coincide, in general, with the standard formulation. 
Morse theory [5] only enables us to say that — under suitable conditions on S[<f>] and O — the 
standard functional integral coincides with an integral over a particular combination of thimbles 
^2 a n a J^<j, for some integer n a , and where the sum runs over all the stationary points of the 
complexified actiorj^J 

However, we have shown in [Tj that the model defined by Eq. ([2]) has the same degrees of 
freedoms, the same symmetries and symmetry representations, and also the same perturbative 
expansion as the standard formulation. In this sense, by universality, the formulation of Eq. ^ 
may be seen as an alternative regularization of the model in Eq. 0. 

3. The Aurora algorithm 

In this section we describe an updated form of the algorithm proposed in [lj. What we want to 
compute are expectation values like: 




where the imaginary part of the action Si is constant in J§ and can be taken out of the integral. 
Hence the phase of e~ s cancels from all observables. Moreover, the real part Sr is bounded 
from below by its value at the stationary poinlrl The functional integral in Eq. Ml) possesses 



2 An argument in [5] (see Sec. 3.2), can be used to show that only the stationary points near to the global minima 
in C can bring a non negligible contribution. See [T], Sec. II. B. 3. 

3 The integral is convergent, because S[4>] is a polynomial in and j7o is union of the curves of SD. 



a bounded real action and hence can be studied, e.g., with a Langevin algorithrrj^J The choice 
of the Langevin algorithm is particularly convenient because the drift term of the Langevin 
equation for the real part of the action coincides with Eq. ([3]), i.e., the SD. Hence, such drift 
preserves the manifold J§ by construction. The Langevin noise, however, needs to be projected 
on the tangent space to Jq, in order to explore the integration domain correctly. 

The projection on the tangent space of J7b at 4> is challenging, because we seem to have 
no way to tell which configurations in the neighborhood of <f> 6 J7q will eventually also fall in 
(/'glob) by following the SD flow. However, the tangent space at the stationary point <?i>giob can 
be computed, by computing the Hessian matrix^ In order to use this fact to generate a noise 
vector tangent to J§ at <f>, we need a way to transport a vector 77 along the flow of SD, while 
keeping it tangent to J$. One way to achieve this is by requesting that the vector 77 is parallel 
transported along the flow of Sr, i.e. requesting that its Lie derivative along dS R is zero: 

CasM = [dS R , V }=0 

which nicely translates into a first order ordinary differential equation, which is also linear in 77: 

J^j( r ) =^2vk{r)dkdjS R . (5) 

k 

Note that we use j, k as multi-indices for (R/I, a, x). 

Actually, we can do better. Eq. @ ensures a parallel transport of 77 along the flow dS R . 
This is actually more than what we need, because we only need to keep 77 tangent to Jq. 
On the other hand, it would be nice to transport 77 by means of an orthogonal flow, that 
preserves orthonormality, isotropy, and also enables the use of stable numerical integrators (see 
Sec. 5.3 of p 7 ]). This can be done quite easily by employing the Iwasawa decomposition of 
8 2 Sr (see, e.g., Sec. VI of [8]), which, at the group level, is a generalization of the Gram- 
Schmidt orthonormalization procedure, and it can be implemented also at the algebra level. 
More precisely, if we write: 

(d 2 S R ) ij = A i>j +D itj + N id , 
where A is skew-symmetric, D is diagonal and N is upper triangular, then the flow 

AiKr) = nir)A{T) (6) 

defines an orthogonal flow whose solutions 77 span the same space spanned by the solutions of 
Eq. ^. This means that the tangent space of J7b is still preserved. But now, the flow defined 
by Eq. ([6]) also preserves orthonormality, isotropy, and such properties can be even ensured 
numerically by using suitable numerical integrators, such as the implicit midpoint rule [TJ. 

We can now summarize the algorithm as follows. We use t to represent the Langevin (Monte 
Carlo) time and r to represent the parameter of the SD. Let us assume that we have already 
a configuration 4>t in j7o, and we also have a path </>t(r) that fulfills the equations of SD and 
connects <j)t to a configuration <j>^ with norm less than e. The value of e must be sufficiently 
small so that the quadratic approximation of the action is valid. One step of Langevin proceeds 
as follows: 

• apply the Langevin force, that consists in taking one step forward along the curve of SD 
which is already available. The result is the configuration (f>' t ; 



3.1 



The phase that comes from the measure d<j) will be treated with reweighting, as discussed in Sec. 
Which is easy to compute also analytically, when g i o b = 0, and can be computed numerically, once and for 
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extract a Gaussian noise rj; 

project out the components of r/ orthogonal to Jq at the origin. More precisely, set 
m = (P — 1)77, where P = d 2 SftJ \J (d 2 Sn), and then re-scale the norm of r/11 according 
to the x distribution; 

evolve 7711 with Eq. (6) from <p[ £ ^ to <p' t . The orthogonal properties of the flow can also be 



preserved numerical 



y by employing, e.g., the implicit midpoint rule. This leads to solve a 



5-D linear system in 7711 (r); 

J t+dt 



define a tentative new configuration as 4*T+dt = + v^tv( T )'i 



use ^t+dt as a starting point to find a solution of the SD equation. This can also be 

integrated with the implicit midpoint rule with the constraint that (P — l)4>f+ dt = 0. This 
leads to a non-linear 5-D system that can be solved — if the initial guess is sufficiently 
good — with a number of Newton-Raphson iterations. 



3.1. The residual phase 

As noted in the previous section, the integral of Eq. Q also includes a residual phase that comes 
from the determinant of the tangent space of Jq at (j>, and is not necessarily real and positive. 
We argue that this should not represent a sign problem. In fact, a sign problem means that the 
average phase \ (Y\ d<j))\ -C 1. However, that phase is constant in the region where the integrand is 
not very small. In fact, such phase is completely neglected in the saddle point expansion, which, 
on the other hand, can be identified with some form of perturbative expansion of the QFT. 
Saying that the average phase could be negligible is equivalent to saying that the "perturbative 
contribution" (in some sense that we cannot define precisely) could be negligible. We do expect 
that a non-perturbative contribution may be important, but that the perturbative one could be 
negligible is contrary to the general expectations. 

In this sense, we have good reasons to hope that residual phase can be treated with 
reweighting, but we have no justification to neglect it. If so, we have to compute it, which 
is quite expensive. This can be done as follow. If we call the complex (2V x 2V) matrix that 
spans the tangent space of Jq a t we have: 

,00 . N T N T N R 

det(T, T ) = - / ds Tr [ -T 0s l^Tr^AT^]^^ AT 0i 

Jt i=0 t=0 k=l 

where the k = 1, . . . , Nr are noisy estimators for the trace. Assuming that the inversion 
of T<j> would require iterations, the computation of the determinant would cost as much as 
evolving NrN vectors, which is described in the previous section. 

Some saving may also come from the fact that the computation of the phase is necessary only 
for the configurations that are actually used for measurement. 



3.2. Preliminary tests 

As a first test, we implemented the method by using the naive Euler integration method. 
Although bound to fail because of the nature of the Euler method, this is an interesting test in 
order to get a feeling of how difficult it is to keep the system on a thimble and how sensitive 
the system is to perturbations. Some results are shown in Fig. [T]and[2] These are encouraging, 
because it seems that even the Euler method is able to keep the system on the thimble long 
enough to see a thermalization. This lets us hope that a better integrator will have good 
chances to converge. Observables appear around the expected order of magnitude, but their 
thermalization is not clear. 



mu=0.9 lambda= 1 .0 mass= 1 .0 

400 r 




Figure 1. The real part of the action on 
a V = 4 4 lattice with Euler integrator. The 
system is stable long enough to show signs of a 
thermalization. 
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Figure 2. Same as Fig. [I] but for a longer time. 
The Euler integrator cannot constrain the system 
on the thimble, and the system drifts away at some 
point. This is also confirmed by the deviation of 
the imaginary part of the action from a constant. 



4. Conclusions 

We have illustrated an new approach to deal with the sign problem that afflicts a class of QFTs. 
It consists in regularizing the QFT on a Lefschetz thimble. Although it does not coincide with 
the usual regularization, it is a legitimate one on the basis of universality. In fact, we could prove 
that QCD on the thimble has the same symmetries, d.o.f. and PT as the usual formulation. 

We have also introduced a Monte Carlo algorithm to achieve an importance sampling of the 
configurations on the thimble. Its numerical implementation will be certainly challenging and 
expensive, but all the steps of the algorithm are, a priori, feasible and have acceptable scaling. 
The residual phase should not give a sign problem (unless we believe that the perturbative 
contribution can be negligible) and hence should be manageable with reweighting, but this 
must be checked. We are presently testing the method for a scalar QFT on tiny lattices with 
encouraging results. 
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